  function ocell = out_moments( G1 , HH , SDX , ZZ , CC , sercell ,shocell,  in ); 
% function ocell = out_moments( G1 , HH , SDX , ZZ , CC , sercell , shocell, in); 
% 
% COMPUTE SIMULATED MOMENTS OF THE SOLUTION FROM GENSYS 
% Code is written for a generic ZZ matrix and cell of matching dimension 
%
% INPUTS
% --------
% Output from Gensys  ZZ is the observation matrix 
% SERCELL   
% IN         Structure with inputs for the simulation 
% 
% in.Ttake    
% in.Tdrop 
% in.nrep        Number of repetitions 
% in.bpind       Index of series to band-pass in the ZZ matrix 
% in.flag_bp  
% in.nk          Number of terms for the filter 
%
% 
% Output 
% -------
% Cell contains std, varf, autocorrelation and correlation of data 
% Showing medians and percentiles 
% 
% Also computes the variance decomposition, although not for filtered data 
% 
% Alejandro Justiniano 9/26/2005 
% ======================================================================
nz = size( ZZ , 1 ); 
nx = size( SDX , 1 ); 

[vdcom,vcfull,sdevf,ACzero,ACone]=vardecom_mode(G1,HH,SDX,ZZ); 

[ee,in] = ch_field( in , 'nrep' , 1000 ); 
[ee,in] = ch_field( in , 'Ttake'   , 150 ); 
[ee,in] = ch_field( in , 'Tdrop'   , 200 ); 
[ee,in] = ch_field( in , 'flag_bp'   , 0); 

Ttake = in.Ttake;
Tdrop = in.Tdrop; 
nrep  = in.nrep; 
flag_bp = in.flag_bp; 

if ee~=0 & in.flag_bp == 1
        
    [ee,in] = ch_field( in , 'bpind' );

    if ee==0 
        error('Must Provide BPIND if wish to compute BAND PASSED series') 
    end 
    bpind = in.bpind;
    
    if length( bpind ) > nz
        error('Bpind exceeds length Z');
    end

    [ee,in] = ch_field( in , 'nk' , round(Ttake^0.25) );
    nk = in.nk;
    Ttake = Ttake + 3*nk + 1;

end 

sercell = sercell(:); 
if length( sercell ) ~= nz 
    error('Mistmatch length SERCELL and zmat matrix')
end 

TT=Ttake+Tdrop;

covmat = zeros( nz , nz , nrep ); 
corrmat = zeros( nz , nz , nrep ); 
aconemat   = zeros( nz , nz , nrep ); 

whandle = waitbar(0,'Computing Simulation'); 

ii=1; 
for ii=1:nrep; 

    waitbar( ii/nrep , whandle ); 
    
    shocks=SDX*randn(nx,TT);

    Ysim=simgensys(shocks,G1,CC,HH,ZZ);
    Ysim=Ysim(:,Tdrop+1:end);
    Ysim=Ysim';
   
    if flag_bp == 1; 
        
        %disp(['Series ',num2str( bpind(:)' ),' will be BP Filtered'] ); 
        %nk = round( size(Ysim,1)^0.25 ) ; 
        Ynew = bpf( Ysim(:,bpind) ,  in.lowp  , in.highp , nk); 
        Ysim( nk+2:end-nk-2 ,bpind)  = Ynew( nk+2:end-nk-2,:);
        Ysim = Ysim( nk+2:end-nk-2 ,:) ; 
        
    end 
     
    covmat(:,:,ii) = cov( Ysim ); 
    corrmat(:,:,ii) = corrcoef( Ysim ); 
    aconemat(:,:,ii) = accov( Ysim , 1);
    
end 
close(whandle); 
T = size( Ysim ); 

% -------------------------------
%           Output cells 
% --------------------------------
pmat = round( nrep*[0.05;0.95] ); 

vcfull = squeeze( median( covmat , 3 ) );
covmat= sort( covmat , 3 ); 
covmat_1 = squeeze( covmat(:,:,pmat(1) ) );
covmat_2 = squeeze( covmat(:,:,pmat(2) ) );

% ==============================
%      Standard Deviation      % 
% ==============================
sdevf  = sqrt( diag( vcfull ) ); 
sdevf_1 = sqrt( diag( covmat_1) ) ; 
sdevf_2 = sqrt( diag( covmat_2) ) ; 
clear covmat*; 

% ------------------------------------------------
%                     titcell: To put at top 
% ------------------------------------------------
ncol = max( 4 , nz + 1 ); 

titcell = emptycell( 5 , nz + 1 );
titcell(1,1) = {'Simulated Moments'};
titcell(1,2:3) = {['N Sim:',num2str(nrep)],['T: ',num2str(T) ]};
titcell(2,1)   = {'Series'}; 
titcell(2,2:nz+1) = sercell'; 

temp_ind = zeros( 1, nz );
if flag_bp == 1 
    temp_ind(bpind) =1 
    titcell(4,1) = {'BP Settings'}; 
    titcell(4,2) = {['NK ',num2str(nk)] }; 
    titcell(4,3) = {['LOWP: ',num2str(in.lowp)] }; 
    titcell(4,4) = {['HIGHP: ',num2str(in.highp)] }; 
end 
titcell(3,1)  =  {'BP Filtered?'};
temp_ans = num2ans(temp_ind);
titcell(3,2:nz+1) = temp_ans';
clear temp_ind temp_ans; 

% =================
% Std and Variance 
% =================
topcell = emptycell(  6 , nz + 1 ); 
topcell(1,1)={'STD'}; 
topcell(1,2:end) = sercell'; 
topcell(2,2:end) = num2cprec( sdevf' );
topcell(3,1) = {'5%'}; 
topcell(3,2:end) = num2cprec( sdevf_1' );
topcell(4,1) = {'95%'}; 
topcell(4,2:end) = num2cprec( sdevf_2' );
clear sdevf_1 sdevf_2; 

topcell(5,1)={'VARIANCE'}; 
topcell(6,2:end) = num2cprec( (diag( vcfull ))' ); 

ocell = merge_cells( titcell, topcell , 1); 
clear titcell topcell; 


vcell = emptycell( nz + 3, nx + 1); 
vcell(1,1:2)={'VARIANCE','DECOMPOSTION'}; 
vcell(2,1) = {'MODEL STA'}; 
vcell(2,2)  ={'Shocks'}; 
vcell(3,2:end) = shocell'; 
vcell(3,1) = {'Series'}; 
vcell(4:end,1) = sercell; 
vcell(4:end,2:end) = num2cprec( vdcom ); 

ocell = merge_cells( ocell, vcell , 1); 
clear vcell; 

% ===============
% Autocorrelation 
% ===============
aconecell = out_sqmat( aconemat, [] , [] , sercell , [] , 'AUTOCORR ONE' ); 
ocell = merge_cells( ocell , aconecell , 1);
clear aconemat; 

corrcell = out_sqmat( corrmat, [] , [] , sercell , [] , 'CORRELATION' ); 
ocell = merge_cells( ocell , corrcell, 1); 
clear corrmat; 




